library(dplyr)
library(stargazer)
library(estimatr)
library(reshape2)
library(texreg)
library(tidyverse)

lag.new <- function(x, n = 1L, along_with){
    index <- match(along_with - n, along_with, incomparable = NA)
    out <- x[index]
    attributes(out) <- attributes(x)
    out
}

lead.new <- function(x, n = 1L, along_with){
    index <- match(along_with + n, along_with, incomparable = NA)
    out <- x[index]
    attributes(out) <- attributes(x)
    out
}

######
## DiD of effect of legalizing same sex marriage on public opinion
######
states2<-readRDS(file="state_covariates.rds")
states2<-states2[,2:4]
states2<-subset(states2, year==2014)
states2<-select(states2, -year)
states2<-rename(states2, abb=POAbrv)
AVA_GayMarriage<-read.csv(file="AVA_GayMarriage.csv")
AVA_GayMarriage<-merge(AVA_GayMarriage, states2, by.x=c("State"),
                       by.y=c("Name"),all.x=T)
AVA_GayMarriage<-AVA_GayMarriage[,1:8]

lm(Favor ~ legal + factor(State) + factor(Year), data = AVA_GayMarriage) %>%
    tidy %>%
    filter(term == "legal")
lm_robust(Favor ~ legal + State + factor(Year), data = AVA_GayMarriage) %>%
    tidy %>%
    filter(term == "legal")
lm_robust(Favor ~ legal + State + factor(Year), data = AVA_GayMarriage,
          clusters = State) %>%
    tidy %>%
    filter(term == "legal")

AVA_GayMarriage$Favor2 <-
    AVA_GayMarriage$Favor/(AVA_GayMarriage$Favor+AVA_GayMarriage$Oppose)

model_estimates<-read.csv("gaymarriage_dyn_psraking8k_dgirtout.csv")
model_estimates<-rename(model_estimates, median=value_mean)
model_estimates2<-read.csv("gaymarriage_disaggregated.csv")
model_estimates$abb<-gsub("abb", "",model_estimates$abb)
model_estimates2$abb<-gsub("abb", "",model_estimates2$abb)
model_estimates<-merge(model_estimates, model_estimates2, by=c("abb", "year"),
                       all.x=T, all.y=T)
AVA_GayMarriage<-rename(AVA_GayMarriage, year=Year)
AVA_GayMarriage<-merge(AVA_GayMarriage,model_estimates,  by.y=c("abb", "year"),
                       by.x=c("abb", "year"), all.x=T)

AVA_GayMarriage<-arrange(AVA_GayMarriage, State, year)
AVA_GayMarriage<-group_by(AVA_GayMarriage, State)
AVA_GayMarriage<-AVA_GayMarriage %>%
    mutate(Favor_lead = lead.new(Favor, 1, along_with = year),
           GM_unweighted_lead = lead.new(naes2004_stategaymarriage_unweighted,
                                         1, along_with = year),
           GM_weighted_lead = lead.new(naes2004_stategaymarriage_weighted, 1,
                                       along_with = year),
           dgo_lead = lead.new(median, 1, along_with = year))	

cor(AVA_GayMarriage$Favor2[AVA_GayMarriage$year==2015],
    AVA_GayMarriage$median[AVA_GayMarriage$year==2015], use="complete.obs")


fit_feedback_model <- function (data, yvar) {
    data$y <- data[[yvar]]
    mod <- data %>%
        mutate(Year = factor(year)) %>% 
        lm_robust(y ~ legal + State + Year, data = ., clusters = State)
    mod$outcome <- yvar
    mod
}

yvars <- c("Favor_lead", "GM_unweighted_lead", "GM_weighted_lead", "dgo_lead")

reg.ls <- plyr::llply(yvars, function (yv) {
    AVA_GayMarriage %>%
        mutate(Favor_lead = Favor_lead / 100) %>%
        filter(year==2014 | year==2015) %>%
        fit_feedback_model(data = ., yvar = yv)
})

texreg(reg.ls, dcolumn = TRUE, include.ci = FALSE,
       custom.coef.map = list("legal" = "SSM Legal"),
       digits = 3, use.packages = FALSE, include.adjrs = FALSE,
       include.rmse = FALSE, table = FALSE)
